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Abstract 


Experimental Gamma-ray attenuation technique along with efficient tomographic 
inversion algorithms and two-fluid modeling, both are in use, to study the liquid metal 
two-phase flow intricacies since long. 

I he formulation lor the chord average void-fractions from experimentally 
collected photon counts (Gamma-ray attenuation technique) and then their inversion to 
have line-by-line, point-by-point and simultaneous void distributions have been 
presented. Seven series expansion methods (Simple ART/Anderson ART/Gilbert 
AR I /Gordon ARI GBH MART/GH M ART/Lent MART) along with two optimization 
reconstruction techniques (MAXHNT/MHM) have been implemented in FORTRAN- 90 . 
The performance of these methods has been checked over three field functions namely 
constant, step, and impulse. Five errors (average error/RMS error/normalized RMS 
error 'maximum absolute error normalized absolute error) have been reported. Plots of 
RMS error vs. relaxation parameter for all the algorithms except maximum entropy 
method have also been presented in step field reconstructions. The reconstructions 
corresponding to the relaxation parameter causing minimum RMS error have been 
presented. 

liven with implications such as effort, cost and time of above discussed 
methodology (Gamma ray attenuation technique), it only gives the void fraction 
distributions at a particular cross-section without any information about fluid velocities, 
pressure etc. 

Additionally, two-fluid modeling of liquid metal two-phase flows for the purpose 
of ADS (Accelerator Driven Sub-critical Reactor) development has been proposed. 
Formulation of fwo-fluid model with proper constitutive laws for ERE (Lead-Bismuth 
Eutectic EArgon system has been done. This is a six-equation model, having two sets of 
field equations, governing the process of conservation of mass, momentum, and energy. 
These six-equations have been solved by an extension of ICE (Implicit, Continuous, 
Hulerian) technique for single phase flows; this technique was named IMF' (Implicit 
Multi-Field) for two-fluid calculations. Formulations for numerical solution have been 
done using finite difference technique with Hulerian cells and interfacial transfer terms 


in 



have been handled linear implicitly (though in actual sense, non-linearity exists). 
Implementation of the model has been carried out in FORTRAN-90 compiler. 

Interfacial drag transport has been formulated from Ishii et al's work, however no entropy 
production at the interfaces and infinitely small interfacial thickness has been assumed to 
deal with interfacial heat transfer. 

In the present work, a simple geometry has been taken into consideration for 
demonstration purpose with the possibility of the application of the two-fluid model in a 
gas-driven ADS target system. The two-fluid model has been solved for bubbly flow in 
undistorted particle regime. Void-fraction patterns along with gas and liquid velocities 
have been reported at four real time levels. 

At last it can be said that all the above studies have been carried out to show the 
possibility of the application of the two techniques (Gamma-ray attenuation technique 
and Two-fluid modeling ) on a gas driven ADS target system to have an insight of two- 
phase flow dynamics there. 
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Chapter 1 


Introduction 


High-density two-phase flows are becoming important in many applications like 
LMMHD, ADS etc. Void fraction is one of the important parameters, which greatly 
influences the performance of these systems because of its direct relation with pressure 
variations. In LMMHD, two-phase flow is developed in the riser tube by introducing gas 
at the bottom to circulate the liquid metal through entire loop. So to design LMMHD 
accurate void fraction patterns need to be studied carefully over the entire length of the 
riser. Two-phase flows are difficult to analyze in wide range of systems. Void 
distributions are governed by complex interactions between two-phases, external 
boundaries and fluid properties. The distribution also depends upon the flow regimes like 
bubbly, churn, slug and annular. For point-by-point, simultaneous and average void 
fraction calculations at a particular cross-section in LMMHD riser leg, Munshi and 
Satyamurthy have already proposed gamma-ray attenuation technique with tomographic 
inversion techniques J 1, 14, 15, and 16J. 

Gas driven ADS as suggested by many researchers, is one more device among utilizing 


liquid metal two-phase flows. The idea to use accelerator driven systems is not new. Over 
the past few decades, however, the motivation for its development has rapidly changed. 
Originally ADS was developed to produce fissile materials, namely 2y) Pu from 1,8 U or 
m U from '''Th in the Material Testing Accelerators project at the Lawrence Livermoore 


Radiation Laboratory at the end of the forties. This project was abandoned in 1952 when 
high grade uranium ores were discovered in the United States. In addition, the subsequent 


idea of exploiting the spallation process to transmute actinides and fission products 
directly turned out to be ineffective. The proton beam currents required (around 30 mA) ! 
were much larger than the most optimistic theoretical designs for an accelerator could j 


achieve. Further diseo\ cries in accelerator technology (also a product of the Strategic 
Defence Initiative - Star Wars program) have started a new wave of interest in using 
ADS. The accelerator generates high energy (around 1 GeV, current 0 and 30 mA) 
charged particles (e.g. protons) which strike a heavy material target. This bombardment 
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Fig- LI Schematic of a typical ADS reactor 

For efficient production of neutrons, the most effective range of incident particles like 
protons is between 1 to 1,5 GeV and the figure of merit for neutron production is defined 
as the number of neutrons emitted per unit amount of energy deposited in the target [10]. 
Spallation target module is one of the important innovative components of the ADS. It 
constitutes the physical interface between the accelerator and the sub critical core, is 


the production of a very intense neutron source, producing spallation neutrons typically 
in the range of 10 1S per second (a process called spallation) [10, 17]. These neutrons 
enter a subcritical core (often called a blanket) where they can be multiplied. Reactor 
with ADS technology operates under sub-critical condition with high level of inherent 
salety, minimum waste production, and can utilize fertile components more efficiently 

[19]. 

Spallation nuclear reactions occur when an energetic projectile (E/A > 100 Mev/nucleon) 
hits a nucleus. I his reaction can be seen as breaking of a nucleus into smaller parts in 
f°rm of nucleons and clusters. 1 he resulting heavy nuclear residues are also highly 
excited and emit a large number of nucleons (basically neutrons) by evaporation process. 
The Schematic of basic spallation process is given in Fig- 1.1. 
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simultaneously subjected to severe thermal-mechanical loads and damage due to high- 
energy heavy particles (proton beam, spallation neutrons, and spallation products). The 
spallation module design should be based on balanced optimization among neutron yield, 
material properties and thermal-hydraulic performance under the required safety and 
reliability conditions [10]. 

One more important part of developing ADS is choosing a suitable target for desired 
performance. Out of various possible options under severe operational conditions, like 
high heat deposition by proton beam, by spallation products and by spallation neutrons 
Satyamurthy et al proposed lead bismuth eutectic (LBE, 45 % lead and 55 % bismuth) to 
utilize at target [10]. It has low melting point 123 °C and also reduces the hazard ol 
polonium gas, forming lead polonide [19]. 

One can have the same liquid acting as reactor coolant, and spallation target (valid for 
fast core). Alternatively, a separated loop for both of them can also be considered. In 
separated loop systems, two schemes are wide spread, window or windowless target 
systems [26-28]. In the window concept spallation heavy metal target and proton beam 
pipe (vacuum environment) are separated by solid window barrier through which the 
proton passes to interact with the target. The alternative concept is windowless, where the 
interface is the free surface of the liquid metal spallation target. Latter is only possible 
with LBH or Lead, because of their low vapour pressures at the operating temperatures 

( < 50()”( ') [ 1U[. 

High heat deposition in the liquid metal target occurs because of the proton beam 
interactions; therefore this heat should be carried away prior to the safety and corrosion 
limits. Various types of ADS schemes have been proposed consisting of fast, thermal, 
and fast thermal cores with different fuel compositions and different coolant schemes 
[19-25]. Mainly, three schemes are popular, natural convection, forced convection 
(mechanical pump driven conditions), and buoyancy driven enhanced by gas injection 
[22]. Above mentioned schemes together with their implementation and operating 
complicacies have been extensively discussed. 

In buoyancy driven systems volumetric expansion coefficient of lead or LBE with high 
boiling temperature helps to circulate the lead or I, BE. The heat deposition by proton 
beam raises the temperature of spallation zone and temperature of cold leg is maintained 
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using a heat exchanger. Higher the temperature difference between spallation zone and 
cold leg, higher is the flow rate, but if temperature difference exceeds a certain limit then 
it causes severe problem of corrosion in the module [29]. Though these systems are 
passively safe but at the same time they have two main disadvantages, namely, required 
target loop size is very large and shooting of coolant temperature at start up. To get rid of 
over shooting at the start up gas injection can be used. Fig-1.2 demonstrates simple 
buoyancy driven ADS target system. 

Coupling flange 

Coolant outlet 
Coolant Inlet 

Proton beam pipe 

Heat exchanger 

Hot rising liquid 
metal 

Cold descending liquid 
metal 



Window 

Spallation region 
Flow guide 


Fig- 1*2 Buoyancy driven target module 

Although, high flow velocities of heavy liquid metals, limits the degree of compactness, 
forced convection systems can be developed with lots of flexibility and compactness. EM 
(electromagnetic) pumps are more favorable than ordinary centrifugal pumps. 

Third scheme is Buoyancy driven enhanced with gas injection systems (Fig-1.3) [27], An 
inert gas is injected to enhance the buoyancy driven flow in these systems. The injected 
gas creates density difference between the two legs that leads to an induced pressure 
difference, which helps to circulate the liquid metal through target module. A heat 
exchanger is used to under limit the operating temperature of cold leg. The inert gas is 
separated out in separator tank located at the top of the module. This system reduces the 
liquid metal inventory and the size of natural circulation system. The need of 
maintenance of the inline dynamic machinery like pump is avoided and also the risk of 
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temperature overshooting at the start up is reduced. Fig- 1.3 gives the details of a gas 
enhanced ADS target module. 


Gas 


.... Coupling flange 

Coolant out 
- Coolant in 
Separator 


;r4 



Proton beam 
Heat 

exchanger 


visai 


Mixer 


Window 

Spallation 

region 

Flow 

guide 


Fig-1.3 Gm enhanced or gas-driven target module 

The studies of liquid metal two phase flows are carried out keeping in mind the 
development of gas enhanced buoyancy driven window type ADS. The performance of 
these systems crucially depends upon the two-phase flow dynamics and void fraction 
distributions. 


Two of various available techniques for two-phase flow studies are considered here, first 
is the C lamina ray attenuation method to calculate chord average void-fractions then their 
inversion using tomographic techniques ( ART/M ART/MENT/MEM) in point-by-point 
and simultaneous void distributions. (I, 30-45 J. This method suffers certain limitations; it 
can give void-fractions at a particular cross-section without any knowledge about 
pressure and velocities moreover experimental projection data is required, which is an 
expensive and time taking activity. 

Additionally, two-fluid model has been used to have an insight of actual liquid metal 
two-phase flow fluctuations {2,8 j. 
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1.1 Review of literature 


A review of the existing literature on studies related to ADS, tomographic techniques, 
and two-fluid modeling has been presented here. First part deals with the on going 
research over ADS technology, secondly advancements in tomographic techniques and 
lastly two-phase flow modeling. 

Satyamurthy and Venkatramani (2002) [10] presented a state-of-art review of the ADS 
technology and developments. ADS studies have been taken up in view of their inherent 
safety and ability to utilize thorium and to transmute present nuclear waste. The main 
thermal-hydraulic feature of the target, heat deposition by proton beam has been 
thoroughly studied. Various issues related to thermal-hydraulics of window and 
windowless configurations based on buoyancy as well as gas-driven systems have been 
presented. The parameters of the target of 1.75 MW (proton beam of 350 MeV and 5 mA 
current) of the demonstration ADS reactor based on quasi-one-dimensional codes have 
been done along with the necessity to develop 2-D/3-D flow models for spallation flow 
region and the free surface flow in the case of windowless configuration. Finally, the 
conceptual design of an experimental LBE facility to be set up for validating the CFD 
codes for target loops for both window and windowless configuration has been presented. 

Satyamurthy and Biswas (2002) |17| presented an interesting study on the thermal 
hydraulics of ADS systems. Preliminary thermal-hydraulic studies carried out on the 
spallation target of window type ADS based on LBH, with a thermal power of 100 MW 
in the fast /one have been discussed. The studies have been carried out both for buoyancy 
as well as gus driven target systems and complete analysis of the loop consisting of 
spallation region, riser, and two-phase flow region of the riser (for gas-driven target) heat 
exchanger, downcomer etc has been done. Parametric analysis for different flow rates, 
beam current (2 and 3 mA and 1 GeV), geometry etc have been presented and 
comparison has been made between buoyancy and gas-driven systems. 

Radon (1917) [50], gave the idea of reconstructed a function from its projections. In the 
strict sense of word, a projection at a given angle is the integral of the image in the 
direction specified by that angle. 
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Hounsfield’s (1973) [39] invention of the x-ray computed tomographic scanner for which 
he received a Nobel Prize in 1972 [4]. Hounsfield used algebraic techniques, like 
ART/MART. 

Ramachandran and Lakshminarayan (1970) [38] developed a new method Convolution 
Back Projection. These algorithms considerably reduced the processing time for 
reconstruction, and the image produced was numerically more accurate. Convolution 
eliminates the step of computing Fourier transform and then its inverse so requires less 
computing time. 

Gordon et al (1970) [33] presented an entirely different approach for tomographic 
imaging consists of assuming that the cross-section consists of an array of unknowns, and 
then setting up algebraic equations for the unknowns, in terms of the measured projection 
data. These techniques are very helpful when it is not possible to measure large number 
of projections and the projections are not uniformly distributed. These techniques are 
known as Series Expansion Methods or Algebraic Reconstruction 
Techniques 'Multiplicative Reconstruction Technique. 

Kak and Slaney in their book “, Principles of Computerized Tomographic Imaging , IEEE 
Press, 1988”, presented a brief discussion over various tomographic techniques. 
Algebraic reconstruction techniques (ART) are discussed in detail in its 7 th chapter and 
suggested, although conceptually this approach is much simpler than the transform based 
methods, for medical applications it lacks the accuracy and speed of implementation. 
However, there are situations where it is not possible to measure a large number of 
projections, or the projections are not uniformly distributed over 180 or 360 techniques 
like ART can be used effectively. 


Gull and Newton (1986) [31] suggested many functions that can be extremized to have a 
realistic solution of reconstruction problem. The two functions that look attractive are 
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entropy and energy. First gives maximally unbiased and latter gives maximally biased 
solution. 

Subbarao at el (1996) [30] thoroughly studied the performance and convergence of 
various series expansion and optimization reconstruction techniques. Algorithms used are 
mainly divided in three groups, namely the additive algebraic reconstruction technique 
(ART), the multiplicative algebraic reconstruction technique (MART) and the 
optimization reconstruction techniques. Additive ART algorithms have shown a 
systematic convergence with respect to number of projections and the value of relaxation 
parameter. MART algorithms have been found producing less error at convergence 
compared to additive ART but converge only at small values of relaxation parameters. 
Increasing noise level in projection data increase the error in reconstructed field and 
decrease the convergence rate. 


Biswas and Munshi (2001) [1] applied ART algorithms over experimentally collected 
projections on I.MMHD experimental facility available in BARC. First ART algorithms 
had been implemented over simulated data thereafter experimentally collected corrupted 
chord average void fractions were inverted to have void distributions. 

Stewart and Wendroff (1984) [46] presented the basic principles of two-phase flow 
modeling and discussed about the ill posed ness of various wide spread models. 

Ishii (1975) [2} presented an interesting formulation of two-fluid model. It was said that 
two-phase flow suffers interfacial discontinuities to over come these he presented 
interfacial jump conditions. Later in 1984, Ishii and Mishima [8] had presented interfacial 
drag formulations and interfacial area calculations with their general applications. 

Liles and Reed (1978) [51] developed a new technique, semi-implicit method, for solving 
the equations of two-phase flow dynamics. Although its applicability over several sets of 
equations representing two-phase flows, including two-fluid model had been proposed, 
but they demonstrated the numerical technique for one-dimensional drift-flux model. Few 
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advantages of this technique over other techniques had also been reported, like stability, 
ease for complicated flow networks, and ease of extension to problems in two or three 
dimensions. 

Harlow and Amsden (1971) [3] revised the ICE (Implicit Continuous-field Eulerian) 
technique and generalized for the arbitrary equation of state and full viscous stress 
tensors. It was said that this method is useful for the numerical solution of time- 
dependent fluid flow problems in several space dimensions, for all Mach numbers from 
zero (incompressible limit) to infinity (hypersonic limit). 

Hirt, Amsden and Cook (1974) [51] given a brief theory of a technique similar to ICE 
with finite difference formulations over three kind of cells, namely Eagrangian (vertices 
can he moved with the fluid), Eulerian (to be held fixed) and be moved in any other 
prescribed manner. Their publication describes the basic methodology, finite difference 
approximations, and discussions about stability, accuracy and zoning. 

Harlow and Amsden (1975) [4] had presented the simplified and extended form of IMF 
(Implicit Multi-field Technique, an extension of ICE for two-phase flow calculations) on 
momentum coupling among fields. The extensions included the capability for 
compressibility in each material phase, the addition of more interpenetrating fields, and 
allowance for the motion of liquid or vapor through a close-packed field of particles. The 
technique had been illustrated by computer generated plots from time varying three-field 
calculation in a c\lindrically. 

Harlow and Amsden (1975) [7] proposed a new technique analogous to ICE for 
numerical solution of two-fluid model (IMF), in which several fields interpenetrate and 
interact with each other. An implicit coupling for each field between mass transport and 
equation of state allowed calculations in all flow-speed regimes, from far subsonic 
(incompressible) to far supersonic. In addition, the momentum transport between fields 
was coupled implicitly, allowing for all degrees of coupling, from very loose to 
completely tied together. Considerable generality has been presented, to permit 
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application to wide scope of complicated problems, for example, the fluidized dust bed, 
the flow of a liquid with entrained bubbles, and atmospheric condensation with the fall of 
precipitation. 

Stewart (1979) [5] presented a study of stability on two-fluid modeling, it was said a two- 
fluid modeling may yield a system of partial differential equations having complex 
characteristics; this results in a mathematically ill-posed initial-value problem. Despite 
that fact that no finite-difference method for solving such a problem can be stable in the 
usual sense, finite difference solution of the two-fluid model is in wide spread use. He 
investigated the numerical behavior of one such set of difference equations, derived 
conditions under which solutions appeared to be well behaved, and offered a physical 
interpretation. 


1.2 Objective of the present work 

Two-phase flow dynamics crucially depend upon the void fraction distributions. Munshi 
et al (1997) [1, 15] suggested iterative tomographic algorithms for the inversion of 
corrupted pi ejection data, collected using gamma ray attenuation technique to calculate 
line-by-line, point-by-point and simultaneous void fraction distributions. A FORTRAN- 
90 code has been developed for series expansion and optimization reconstruction 
algorithms, to check their performance over three simulated functions namely constant, 
step, and impulse with a motive to apply these algorithms over ADS leg. Mainly five 
errors have been calculated in order to check their performance for the application over 
an ADS leg, the range of relaxation parameter of all the algorithms with step 
reconstructions has been reported, different relaxation parameter reconstruct different 
features of a simulated data. If one algorithm shows divergence for some value of 
relaxation parameter then it can be concluded that the algorithm is not capable to 
reconstruct the features corresponding to the concerned relaxation parameter. The 
features corresponding to minimum RMS error are most desirable for better locally 
spread reconstructions thus all the reconstructions have been presented corresponding to 
the relaxation parameter having minimum RMS error. 
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For the same two-phase flow studies, two-fluid model has been formulated to study the 
liquid metal two-phase flow’s complex microphysics in pure bubbly flow under viscous 
regime with appropriate analytical models for drag force [2, 8, 12]. The finite difference 
discretized equations have been solved using IMF (implicit multi-field technique) 
technique, an extension of ICE (implicit continuous technique) of single-phase flows [3, 
4]. Convection terms have been discretized using partial donor cell differencing. A 
FORTRAN 90 code is developed for the numerical solution. 

In the present work six field equations are solved for demonstration purpose over a 
simple geometry in bubbly flow under viscous regime. Liquid and gas velocity vectors 
along with void fraction contour are plotted. 

The tomographic algorithms and transient two-dimensional two-fluid model have been 
developed in the present work, keeping in mind the two-phase llow studies for the design 
and development of gas-driven accelerator driven systems. 

1.3 Structure of the thesis 

Chapter 2 presents an analytical study of projection data calculations from experimentally 
collected photon counts and discussions about various tomographic algorithms. Chapter 3 
deals with the mathematical formulations of the possible algorithms for incomplete data 
reconstructions (series expansion and optimization algorithms). In Chapter 4 the 
reconstructions of three simulated projection data with five errors, profiles of RMS error 
vs. relaxation parameter, and discussions over results have been carried out. Chapter 5 
has been devoted to the formulations of transient two-dimensional two-fluid model along 
with interfacial drag, discussions about stability, and applicability on a LBE-Argon 
system. Chapter 6 has been devoted to results and discussions over two-fluid simulations 
of an LBE-Argon system; velocity vectors and void fraction profiles have also been 
presented. At last, in Chapter 7 conclusions and future scope of the present work have 
been presented. 
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Chapter 2 


Reconstruction algorithms and theoretical approach 

Image reconstruction from projection can be viewed as linear inverse problem with 
discrete data. Reconstruction can be achieved mainly by four ways, direct matrix 
inversion, transform methods (Fourier and Hilbert transform), filtered convolution back 
projection methods, and ART. Among all discussed Convolution Back Projection 
Algorithm is used extensively [37-38]. Convolution Back Projection requires data for a 
large number of angles and should be uniformly distributed to have an exact inversion. 
However in many practical cases it might not be possible to have such large data. On the 
other side, with limited data tomography encounters ill-posed problem in reconstruction. 
To overcome these tomographic problems Series Expansion Methods have been 
proposed, like Additive ART/ Multiplicative ART [33, 39 and 40]. As in the conventional 
ART methods, additive or multiplicative corrections are possible in all the algorithms, 
which belong into ART family. Recent studies have indicated that, of the several classes 
of reconstruction algorithms applicable to limited data, those based on multiplicative 
algebraic reconstruction (MART) are the fastest, most flexible and accurate. Verhoeven 
presented various types of MART algorithms and applied them for interferometer data 
[ 40 ]. 

Reconstruction may be defined ill posed based upon two logics, incomplete data and 
noisy data. Being ill-posed problem in reconstruction with incomplete and noisy data, 
there are many feasible solutions, which are consistent with the data. It is a matter of 
practical necessity to select a single solution from the feasible set. The Maximum 
Entropy method consists of choosing a single feasible solution, which has maximum 
configurationally entropy, means it is maximally noncommittal about unmeasured 
parameters. Gull and Newton suggested many functions those can be extremized to have 
a realistic solution to reconstruction problem [31]. The two functions those look attractive 
are entropy and energy. First gives maximally unbiased and latter gives maximally biased 
solution. 
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Presently focus has been kept on series expansion and optimization methods. AH the 
methods are iterative in nature, Gordon first proposed ART algorithm with name Gordon 
ART for CT images [33]. Later Gilbert developed a form of ART called simultaneous 
iterative reconstruction technique (SIRT) [34], After this Anderson and Kak combined 
ART and SIRT and gave a new name to this, simultaneous algebraic reconstruction 
technique (SART) [35]. In later 1 990’s Mayinger opted a simplest reconstruction 
algorithms for reconstructions, Simple ART [32]. 

With these additive correction algorithms, multiplicative correction algorithms are also 
proposed for ill-posed problems. First Gordon et al (1970) proposed GBH MART [33]. 
Later Gordon and Herman (1974) proposed a slight difference in correction and named 
this algorithm GH MART [41 J. Lent (1977) proposed another multiplicative algorithm 
Lent MART [36]. And it was reported that multiplicative algorithms are better and 
efficient for ill-posed problems. 


2.1 Theoretical approach for chord average void-fractions 

Projection of a field in a given direction can be seen as the line integral of the attenuation 
coefficient (p) of the two-phase flow field. Interaction of matter with alpha and beta 
particles is a multiple collision process, however with Gamma rays it is only one-shot 
process. This directly leads to the exponential absorption law for Gamma rays, namely. 


/{.V) = 


( 2 . 1 ) 


Here, I(X) is the intensity at thickness X, and p, the total absorption coefficient. J () , is 
initial intensity of the radiation at the entrance into the investigated volume. 
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Fig-2.1 shows the photon count collection at a cross-section on a pipe under liquid metal 
two-phase flow. 
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Fig-2. 1 Projection data collection 


Neglecting the attenuation coefficient of air, / A/ , is the number of counts per unit second 
collected at the detector end when only liquid metal is flowing through the pipe./ 2tI1 , is 


under the presence of mercury-nitrogen two-phase flow in the pipe. /,, , is the counts per 

second when nitrogen or air is only present in the pipe. Whereas gi, Xj, p. 2 , X 2 , j. 13 , X 3 , fi 
. 1 , X.i, are attenuation coefficient and thickness for lead block 1, liquid metal, lead block 2, 
and pipe material respectively. 


Now, 


A' 


= X t + X : + X } + A v 

hn m ~ I Mu — I An 


( 2 . 2 ) 


I m ( A ) — I M,, e 


Pf X( p>.\ y p t X r fh «V / 


(2.3) 


1 ( $ ) = 1 


-///A )p : X : P<X t //, V f 


(2.4) 


I a( ^ = I An ^ 
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(2.5) 
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Now by equations (2.2), (2.3), and (2.4), 


Iu(X) _ 
h<p(X) 


Now by equations (2.2), (2.3), and (2.5), 


h± (_ X ) _ 
I A (X ) 


( 2 . 6 ) 


(2.7) 


Using equations (2.6) and (2.7), 


In 


a r =■ 


In 


f ho{X)' 

7 a / (aq 
f I a{ X) s 
h, (AO. 


(2.8) 


Where, a v . is the chord average void fraction. 

Projection input for the series expansion optimization algorithms can be calculated by 
multiplying the chord average void fraction by chord length. 


1\ = a, = a r x L, 


(2.9) 


Two-dimensional Held is divided into square cells, and cells are numbered in a regular 
fashion from bottom to top. The void fraction is assumed to be constant throughout a cell 
(i.e.f/ j , is the void fraction of j th cell). In the present analysis, source, detector, and rays 

are considered ideal The length of intersection of i lh ray and j th cell, denoted by w u for 
i~l, 2 ...M and j 1, 2 ...N represent the contribution of j 1,1 cell to the total void fraction 
along the i Ih ray. 

Schematic of the discretization of the physical domain has been presented in Fig-2.2, 
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Fig-2.2 Discretization of field 

By accounting the contribution of each cell into the average void fraction along a ray, the 
result is the system of linear equations. 


( 2 . 10 ) 

Or. ' 1 

( 2 . 11 ) 

If M and N were small, conventional matrix theory methods would have been used to 
invert the system of equations. However, in practice M and N are large enough which 
precludes any possibility of direct matrix inversion. Of course, when noise is present in 
the projection data, and when M<N, even for small N it's quite difficult to use direct 
inversion method; some least square methods may have been used. But when M and N 
are large, such methods are also impractical. Some iterative methods like ART/MART 
and optimization methods can be used for the inversion of the above problem. 
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Chapter 3 


Mathematical formulations for reconstructions 

The incomplete projection data leads to the application of series expansion and 
optimization methods for the present work. Additive as well as multiplicative correction 
algorithms have been studied. Optimization methods (maximum entropy and minimum 
energy methods) have also been used and results thus obtained are compared to report the 
best algorithm among all the available algorithms for the constant, step and impulse field 
reconstructions. The MART algorithms have shown their convergence to optimization 
algorithms with a necessary condition that the feasible set of solution, which satisfies the 
projection data, should be nonempty, means projections should not be corrupted [43, 44]. 
Consider a Cartesian grid of rectangular elements, introduced in a filled layer over a 
horizontal plane. Field parameters to be constructed are denoted by./] . j = 1, , N , 

The projection, which is a path integral, is approximated as 

N 

P, = a, = // / = 1 M (3.1) 

/-/ 


The weight function ny, is defined as the length of intersection of the ray with the 
/'" pixel, it can be seen that there are N unknowns that need to be solved from M 
equations. The resulting matrix [ M »,J is highly sparse and rectangular with sizeM xN . 
And the problem can be viewed as developing generalized inverse of the matrix [w,.,], 
while obtaining the field values J) . 
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The statements for algorithms with additive corrections are given first, followed by 
multiplicative correction algorithms and at last optimization methods have been 
presented. 


3.1 Additive ART 

3.1.1 Simple ART 

Simple ART step by step is as follows [32], 
1 . initialization /” e R" is arbitrary. 


do for each iteration k. 
do for each angle of projection <9 . 
do for each ray i g 

a. compute the correction ((a - (aA ) 

b. compute the total value of weight function along each ray. 



i i 


, , , 0? /(l)i'v/> (tZ itljiln-n 

c. calculate the average weight 1 unction 

w ,o 


end do 

do for each ray i 0 

a. update the each field parameter as 

^ifK/ io)cxp (d ) 


f™ = f t M + k 


w, 


ll) 


end do 

do lor each ray i t/ 

a. update the approximate projections using 

A' 

(of ,0)11,0 ~ Wnt.i J 1 

i*i 
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end do 


end do 


a. stop when abs 


//' , 


X 100 < 0.01 


end da 


3.1.2 Gordon ART 


In this method corrections are applied to all the cells through which i‘ h ray 
passes, before calculating the correction for the next ray. Hence, the number 
of rays per angle of irradiation is not important [33], 

1. initialization /° e R" is arbitrary. 
do for each iteration k 
do for each ray i 

a. compute tire correction - {a ) 

b. compute the correction coefficient 

v 

E i 

H',,; 

/"! 

e. apply a correction to each cell / of the test field through which the 
present ray passes as 

/• II Id - ^'i,i ((^<)«/> (d i ) 

j, ~J> +A — — — : 


Old do 


a. stop when abs 






• K 


x 100 <0.01 


h. compute projection numerically 

.v 

(fif l)lht,i ~ ; = f idl 


end do 
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3.1.3 Gilbert ART 


This method is also known as SIRT (simultaneous iterative reconstruction 
technique). In SIRT, the elements of the field function are modified after all 
the correction values corresponding to individual rays have been calculated 

[ 341 . 

1. initialization /“ e R" is arbitrary. 
do for each iteration k 
do for each ray / 

a. compute the correction ((aXx? ~ (tf, ) 

b. compute the correction coefficient 


N 



r i 


end do 

do for each ray i 

a. identify all the rays, A/„ passing through a given cell and the 

corresponding u-,., and ((«,),, A/ , -(a , ),/,,„)• 

b. apply a correction to each cell j of the test field through which the 


present ray passes as 
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end do 
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stop when abs 
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X 100 <0.01 


b. compute projection numerically 





r-l 


i = 1,. 


.J4 


end do 



3.1.4 Anderson ART 

This algorithm combines the ART and SIRT algorithms. Method of applying 
correction is similar to SIRT but the structure is similar to ART [35]. 

1. initialization /° e R" is arbitrary. 
do for each iteration k. 
do for each angle of projection# . 
do for each ray i„ 

a. Compute the correction (( a , 0 \, xp - {a, oh,*, ) 

b. Compute the total value of weight function along each ray, 

.V 

c i ii = yi Wit>. i 
/■I 

end do 

do for each ray /„ 

a. update the each field parameter as 

j. («•» j- "hi t id 

h = // +A 

L ll) 

end do 

do lor each ray i 0 

a. update the approximate projections using 
,v 

(or ~ 7. WiO.I J / h ~ M o 

nl 


end do 
end do 


Iwlt * 
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stop when abs 
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x 100 < 0.01 


end do 
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3.2 Multiplicative ART 

When the correction is multiplicative, the ART is called multiplicative ART 
(MART) [33, 40 and 41], 

1. Initialization f° e R" is arbitrary. 
do for each iteration k 
do for each ray i 

a. Identify all the rays, M CJ passing through a given cell and the 
corresponding w , ; , (a) exp , and (a ■ 

b. Apply a correction to each cell j of the test field through which 
the present ray passes as 


3.2.1 GBH MART 
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3.2.2 GH MART 
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3.3 Optimization techniques 
3.3.1 Maximum entropy method 

The probability of finding a system in a given state depends upon the multiplicity of that 
state. 1 hat is to say, it is proportional to the number of ways you can produce that state. 
Here a "state" is defined by some measurable property, which would allow you to 
distinguish it from other states. A greater value of multiplicity of state implies greater 
possibility of the system in that state. Based on the above-mentioned definition of entropy 
a tomographic algorithm is developed by Gull and Newton [31]. This method produces 
an unbiased solution and is maximally non-committal about unmeasured parameters. A 
formulation for the present algorithm is given as. 

Maximum entropy functional is defined as 




r-l 


(3.2) 


Subject to the constraints 


a, 2_, »'(,(/ i i = h , M (projection data) 

i i 

f t > 0 (a prior condition) 


(3.3) 


The 1 ay.anyian multiplier technique has been used for maximization of the functional F 
under imposed constraints. The problem now reduced to the solution of a set of non- 
linear equations with unknown I ugrangian multipliers. Non-linear equations are solved 
after linearization using Taylor series expansion method, by Gauss-Seidel iterative 
technique (APPENDIX A) [45] 
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3.3.2 Minimum energy method 

When a system is having minimum energy, then system is in most stable state. Based on 
the above definition formulations for minimum energy method is done as follows [31], 


Minimum energy functional is defined as 


// 

/=/ 


Subject to the constraints 


X 

a, = f , i = 1 M (projection data) 

i i 


(3-4) 


(3.5) 


The Lagrangian multiplier technique has been used for minimization of the functional F 
under imposed constraints. Problem reduced to the solution of a set of linear equations 
with Lagrangian multipliers unknown. Standard Gauss-Seidel method is used to solve the 
linear equations in I agrangian multipliers and then filed values have been calculated 
(APPHND1X A) [45 J. 
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Chapter 4 


Simulated Results and discussions on tomographic methods 

The performance of the aforesaid discussed algorithms has been checked over simulated 
data. Mainly three field functions have been used to test performance constant, step, and 
impulse inside a unit circle. 61 rays and 12 views in 180° total view angle have generated 
projection data. Ray spacing during reconstructions has been taken 0.01639. During all 
reconstructions the initial values in the pixels inside the circle are taken 1.0 and outside 
forced to zero for the purpose of the application of the above discussed algorithms on a 
vertical pipe consisting of two-phase liquid metal flows. 


4.1 Reconstruction errors 

During the present work, five error norms have been reported. 
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Where, i s assumed field for simulated projections, and is 

reconstructed field using simulated projections. 

Relaxation parameter in algebraic reconstruction techniques is analogous to the filter 
functions in transform methods [42]. Opting for an appropriate relaxation parameter for 
minimum error reconstructions, in other way the most suitable relaxation parameter to 
watch out maximum features, is very' important. Relaxation parameter has to do nothing 
with the convergence and higher relaxation parameter may take longer time than small 
one. It was assumed that if RMS error is at its minimum in an image reconstruction then 
the reconstructed image is with maximum possible features with that particular algorithm. 
E , (RMS error) is used to find out the relaxation parameter (/l), relaxation parameter 
thus selected is corresponding to the minimum E A (RMS error). Fig-(4.1)-(4.8) shows the 
variation of RMS error with relaxation parameter in step field reconstructions. 
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Fig-4. 1 RMS error vs. relaxation parameter Fig-4.2 RMS error vs. relaxation parameter 

for Simple ART for step field reconstruction for Anderson ART for step field reconstruction 
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Fig-4.3 RMS error vs, relaxation parameter Fig-4.4 RMS error vs. relaxation parameter 

for Gilbert AR1 for step field reconstruction for Gordon ART for step field reconstruction 
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Fig-4.5 RMS error vs. relaxation parameter Fig-4.6 RMS error vs. relaxation parameter 

for G'BIf MARI for step field reconstruction for GH MART for step field reconstruction 
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Fig-4,7 RMS error vs. relaxation parameter Fig-4.8 RMS error vs. relaxation parameter 
for Lent MART for step field reconstruction for minimum energy for step field reconstruction 

Fjg-4. i clearly shows that at the lower relaxation parameter in simple ART, RMS error is 
having steep change, later nearly a linear decrease. This trend also suggests that simple 
ART favors maximum and fastest features in an image reconstruction. Fig-4.2 does not 
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have any abrupt change in values, it rather gives nearly smooth change, and Anderson 
ART can reconstruct limited and slow features. Fig-4.3-4.4 again shows a smooth 
decrease, and suffers rather smaller range of relaxation parameter. Gilbert and Gordon 
ART cannot reconstruct the image with fast changing features, though it is desirable for 
better reconstruction. Fig-4.5-4.6 gives the profile of RMS error with relaxation 
parameter for GBH and GH MART respectively, having multiplicative corrections, and 
can reconstruct limited and slow features because higher relaxation parameter is not 
favorable. Fig-4.7 demonstrates the performance of Lent MART, preferring slow changes 
as of G1 1 MART; error level is also small in this algorithm and most preferred among all 
MART algorithms. Fig-4.8 shows RMS error variation with relaxation parameter, and 
minimum energy method is found producing fast and maximum possible features in 
reconstruction and error level in this algorithm is also at minimum. The field values in the 
present step reconstruction, does not have fluctuating changes, and better image can be 
reconstructed with fast changing features (i.e. highest possible relaxation parameter). 
Table-4. 1-4.21 demonstrates five errors in simulated constant, step and impulse field 
reconstructions. 
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Table-4.1 Different error levels, number of iterations, and relaxation parameter are given for Simple ART 

constant function reconstructions 
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Constant field 
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Table-4,2 Different error levels, number of iterations, and relaxation parameter are given for Anderson ART 

constant function reconstructions 


Constant field 

„ _ 

2 

p 

a 

p 

Ec 

E 0 

E, ; 

Iterations 

0 Dl 

I.225E-02 

4.729E-02 

4.432E-02 

9.712E-03 

1.942E-02 

195 

01 

i 186H-02 

4.578E-02 

4J98E-02 

9.468E-03 

1.893E-02 

55 

1 ~0.l5 

_______ 

4J85E-02 

4,1 12E-02 

8.663E-03 

L732E-02 

126 


Table-4,3 Different error levels, number of iterations, and relaxation parameter are given for Gilbert and 

Gordon ART constant function reconstructions 


Constant field 


A 

k A 

E a 

E c 

f 

e e 

Iterations 

0 01 

2.967E-02 

LI45E-0! 

8.506E-O2 

2.363E-02 

4.726E-02 

355 

u o’ 

I61E-02 

1 219E-0I 

9.17IE-02 

2.502E-02 

5.004E-02 

530 


Table-4,4 Different error levels, number of Iterations, and relaxation parameter are given for GBH MART 

constant function reconstructions 
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Constant field 

A 

Ea 

Eb 

Ec 

e d 

E, 

Iterations 

0.01 

8.677E-03 

3.348E-02 

3.080E-02 

6.942E-03 

1.388E-02 

264 

0,1 

8.527E-03 

3.290E-02 

3.199E-02 

6.788E-03 

1.357E-02 

58 

0 2 

9. i 00E-03 

3.51 IE-02 

3.451E-02 

7.240E-03 

1.448E-02 

49 

0,22 

| 9.934E-03 

3.833E-02 

3.788E-02 

7.902E-03 

1.580E-02 

73 

0.23 

i 110E-03 

4.285E-02 

4.287E-02 

8.831E-03 

1.766E-02 

196 


Tabfe-4.5 Different error levels, number of iterations, and relaxation parameter are given for GH MART 

constant function reconstructions 


r ' - — — — — — — — 

j Constant field 

A 

E, 

E H 

Ec 

e d 

E, : 

Iterations 

001 

8,66011-03 

3.341E-02 

3.075E-02 

6.929E-03 

1.385E-02 

258 

6 i 

K3 IDE-03 

3.206E-02 

3.1 J5E-02 

6.616E-03 

1.323E-02 

57 

0 2 

8 1 73E-03 

3.153E-02 

3.092E-02 

6.51 IE-03 

1.302E-02 

44 

~ 0 22 


| 3 ,0378*02 

3.009F.-02 

6.275E-03 

1.255E-02 

76 

023 


2,67 IE-02 

2.723E-02 

5.528E-03 

1.I05E-02 

212 


Table-4,6 Different error levels, number of iterations, and relaxation parameter are given for Lent MART 

constant function reconstructions 


Constant field 


r' 

F 

e b 

E c 

Ed 

E, 

Iterations 

$1™ 

i 7951 -02 



6.928E-02 

7.785E-02 

I.425E-02 

2.850E-02 

- 

61 

' cTi 


3,1188-02 





0 5 

1 2771 **02 

4 9308-02 


1.015E-02 

2.030E-02 

174 

US 

i 087E-02 

4. 198E-02 

4.049E-02 

8.540E-03 

1 ,7098-02 


0? 

4 172E-03 

1,6108-02 

E545E-02 

3.397E-03 

6. 795 £-03 

963 

O 

"4 097K-03 

I.581E-02 

2.570E-02 

3.226E-03 

6.453E-03 

1459 

I S 

,„„„T 

L508E-02 

4.03 IE-02 


5.880E-03 

2017 

l 8 “ 

[ "4 9028-03 

i 

1.89 IE-02 

8.870E-02 

2.439E-03 

4.879E-03 

7659 


Table-4.7 Different error levels, number of iterations, and relaxation parameter are given for minimum energy 

constant function reconstructions 
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Step field 

x 

Ea 

e b 

Ec 

e d 

e e 

Iterations 

0.01 

6.98 IE-02 

1.755E-01 

2.017E-01 

4.976E-02 

6.863E-02 

206 

0.1 

3.098E-02 

7.79 IE-02 

I.382E-01 

1.937E-02 

2.672E-02 

183 

0.5 

2.995E-G2 

7.534E-02 

1.377E-01 



79 

0.9 

2.985E-02 

7.509E-02 

1.377E-01 



54 

1 3 

2.978E-02 

7.491E-02 

1.376E-01 

1.947E-02 

2.686E-02 

44 

1 8 

2.943E-02 

7.403E-02 

1.346E-01 

1.929E-02 

2.661E-02 

97 

1.9 

2.938E-02 

7.389E-02 

1.341E-01 

1.926E-02 

2.657E-02 

104 


Table-4 J Different error levels, number of iterations, and relaxation parameter are given for Simple ART step 

function reconstructions 




Step field 




A 

E., 

E H 

E c 

E n 

e e 

Iterations 

o' 01 

3 12SE-02 

7.860E-02 

1.405E-01 

1.944E-02 

2.681E-02 

197 

0 I 


7509E-02 

1 376E-01 

1.948E-02 

2.687E-02 

57 

0 5 


7, 206 E -02 

I 228E-01 

1.891E-02 

2.609E-02 

205 


T 82 5 E -02 

7.I05E-02 

1 154E-01 

1 858E-02 

2.562E-02 

254 

fl 

“822H-02 

7 098 E “02 

1. 1 5 IE-01 

1 .859E-02 

2.564E-02 

294 


Table-4*<> Different error levels* number of iterations, and relaxation parameter are given for Anderson ART 

step function reconstructions 


Step field 


I ‘ ' 

E < 

h B 

E c 

E n 

E,.; 

Iterations 

001 

3.0I5E-02 

7.583E-02 

L382E-01 



354 

o.t 

2.987E-Q2 

7.513E-02 

i,362E-01 

1.959E-02 

2.7O2E-02 

82 

0 15 

2 965E-02 

7.458E-02 

L342E-01 

E946E-02 

2.684E-02 

118 


Table-4. IQ Different error levels, number of iterations, and relaxation parameter are given for Gilbert and 

Gordon ART step function reconstructions 
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Step field 

A 

£, 

Eb 

E c 

i 

p 

Ee 

Iterations 

0,01 

4.347E-02 

1.093E-01 

2.000E-01 

3.341E-02 

4.608E-02 

413 

0.02 

4.396E-02 

I.105E-01 

2.032E-0I 

3.388E-02 

4.673E-02 

291 


Table-4. II Different error levels, number of iterations, and relaxation parameter are given for GBH MART step 

function reconstructions 


Step field 


A 

E 

U A 


E c 

E n 

E f . 

Iterations 

0 01 

3 136E-02 

7.888E-02 

1.407E-01 

2.025E-02 

2.793E-02 

505 

o i 

3 105E-02 

7,8 IDE -02 

I.405E-01 

2.061E-02 

2.843E-02 

109 

1 2 

3. i I2E-02 

7 827E-02 

1.408E-0I 

2.070E-02 j 

2.865E-02 

68 

Wn~ ~ 

3.I26E-02 

I.863E-02 

I.4I6E-01 

2.097E-02 

2.892E-02 

70 

' oU” 

! 3 I37E-02 

1 

1.890E-02 

1.408E-01 

2.123E-02 

2.928E-02 

188 


Table-4 J 2 Different error levels, number of iterations, and relaxation parameter are given for GH MART step 

function reconstructions 


Step field 


A 



E c 

U l) 

E i; 

Iterations 

ooi 

3 1 J6H-02 

7 888E-02 

1 405E-01 

2.024E-02 

2.792E-02 

503 

ol 

I 102E-02 

7 802H-02 

I 403H-01 

2.038E-02 

2.839E-02 

109 

0 22 

3 096E-02 

7 787E-02 

1,3981’ -01 

2.06 IE-02 

2.843E-02 

72 

02} 

I 066E-02 

1 71 IE-02 

I.366E-0I 

2.Q45E-02 

2.820E-02 

198 


Table-4, 0 Different error levels, number of iterations, and relaxation parameter are given for Lent MART step 

function reconstructions 
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te p field 


A 

E., 


E c 

E D 

Es 

Iterations 

0,01 

5.587E-02 

1.405E-01 

1.731E-01 

4.351 E-02 

6.001 E-02 

43 

d. i 

3.021 E-02 

7.599E-02 

I.306E-01 

2.064E-02 

2.847E-02 

1289 

6 5 

2.822E-02 

7.098E-02 

U46E-01 

1,861 E-02 

2.567E-02 

2361 

oi 

2.817E-02 

7. 08 5 E-02 

U47E-01 

I.863E-02 

2.569E-02 

2407 

1,0 

2.817E-02 

7.086E-02 

1.150E-01 

1 . 865E-02 

2.573E-02 

2440 

1 3 

2.820E-02 | 

7.090E-02T 

I.152E-01 

1.868E-02 

2.577E-02 

2929 

1.5 

— 

2. 824 E-02 

7,104E-02~ 

1.159E-01 

1.873E-02 

2.583E-02 

4106 


2.850E-02 

7.I69E-02 ~ 

1.176E-01 

1.905E-02 

2.620E-02 

5044 

19 

2. 863 E-02 

7.201 E-02 ' 

1.209E-01 

1.898E-02 

2.618E-02 

12008 


I abIc-4.14 Different error levels, number of iterations, and relaxation parameter are given for minimum energy 

step function reconstructions 



^ m Pul$e field 

A 


E/i 

E c 

E d 

E, 

Iterations 

o.o i 

1,3490*01 

1 .0380-0 1 

6.900E-01 

1.163E-01 

4.646E-02 

167 

o i 

5 j 05 E-02 

3.9300-02 

3.952E-0I 

3.619E-02 

1 .445 E-02 

97 

ps 

4.5620-02 

3,5120-02 

3.022E-0I 

3.3 5 50-02 

1,3 39 E-02 

46 

(19 

4 2581*4)2 

32780-02 

2 674E-01 

3. 1400-02 

1 ,253 E-02 

125 

1 3 

4 008E-02 

3.0850-02 

2 662E-0J 

2,9050-02 

1.I60E-02 

181 

i s i 

3.7920-02 • 

2.9190-02 " 

2.647E-01 

2.71 7E-02 

1.084 E-02 

220 

i'9 “j 

3 763 E-02 

2,8970-02 " 

2.64SE-0 1 

2.6920-02 

1.07 5 E-02 

222 


Table-*.! 5 Different error levels, number of iterations, and relaxation parameter arc given for Simple ART 

impulse function reconstructions. 

















Impulse field 

A 

Ea 

e b 

E c 

e d 

e e 

Iterations 

0.0 1 

4.644H-02 

3.575E-02 

2.950E-01 

3.456E-02 

1.380E-02 

308 

0,1 

4.SS0E-02 

3. 5 03 E-02 

2.717E-01 

3.403E-02 

1.359E-02 

65 

0 15 

4.357E-02 

3.354E-02 

2.684E-01 

i 

3.229E-02 

1.289E-02 

114 


Table-4. 1 7 Different error levels, number of iterations, and relaxation parameter are given for Gilbert and 

Gordon ART impulse function reconstructions 


Impulse field 


1 

e a 

E, 

E c 

&D 

E s 

iterations 

“00! 

2 038E-0I 

I 569E-01 

6.616E-01 

1.605E-01 

6.408 E-02 

392 

' 0 02 

| 2 234E-01 

1 700E-01 

7.200E-01 

1.850E-01 

6.704E-02 

560 


Table -4 .18 Different error levels, number of iterations, and relaxation parameter are given for GBH MART 

impulse function reconstructions 


Impulse field 

" 1 A 

E, 

E„ 

E c 

E» 

E s 

Iterations 

0 01 

6 205E-02 

4.770E-02 

3.353E-01 

4. 794 E-02 

1.914E-02 

404 

FT 

T448E-02 

4, 193 E-02 

2.897E-01 

4.198E-02 

1.6 7 6E-02 

84 

02 

l64ii::o2 r 

2.803E-02 

2.220E-0I 

2.652E-02 

1.059E-02 

48 

0 22 

3.40 IE-02 

2.6I8E-02 

2.037E-0I 

2.433 E-02 

9.715E-03 

52 

0 23 ' 

y.inBii 

2.S34E-02 

2.102E-01 

2.328E-02 

9.297E-03 

115 


Table-4.19 Different error levels, number ofiterations, and relaxation parameter are given for GH MART 

impulse function reconstructions 
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Impulse field 

A 

e a 

Eh 

E c 

E d 

bn 

bn 

Iterations 

0.0 1 

6.227E-02 

4.794E-02 

3.361E-01 

4.813E-02 

1.921E-02 

415 

0,1 

5.105E-02 

3.930E-02 

3.952E-01 

3.619E-02 

1.445E-02 

97 

0.2 

S.908E-G2 

4.548E-02 

3.030E-01 

4.570E-02 

1.827E-02 

78 

0,22 

5.843E-02 

4.498E-02 

3.026E-01 

4.519E-02 

1.804E-02 

85 

0.23 

5.277E-02 

4.Q60E-02 

2.993E-01 

4.025E-02 

1 .607E-02 

232 


Table-4.20 Different error levels, number of iterations, and relaxation parameter are given for Lent MART 

impulse function reconstructions 


Impulse field 


A 

E 

1 A 

E {l 

E c 

e d 

Eh 

Iterations 

« oi 

I 290E-01 

9.936E-02 

6.958E-01 

1.026E-01 

4.099E-02 

27 

0 I 

4.547E-02 

3.500E-02 

2.865E-01 

3.212E-02 

1.367E-02 

1540 

0 5 

3.2 15 1:4)2 

2.475E-Q2 

2.432E-01 

2.201E-02 

8.789E-03 

1040 

0 8 

2 750E-02 

2.1 17E-02 

2.160E-01 

1.680E-02 

6.71 IE-03 

2074 

nr 

2.810E-02 

2.164E-02 

2.208E-01 

1. 7280-02 

6.901E-03 

798 

i 1 


2.120E-02 

2.176E-0I 

1.682E-02 

6.716E-03 

3714 

o 

’ 2.77IE-02 

2.I33E-02 

2.175E-01 

1.692E-02 

6.755E-03 

4275 

1 8 


2.174E-02 

2.203E-01 

1.735E-02 

6.928E-03 

9095 

1 5 

TOToT 

2.184E-02 

2.212E-01 

1.780E-02 

[ 7.109E-03 

14697 


Table-4. 2 1 Different error levels, number of iterations, and relaxation parameter are given for Minimum energy 

impulse function reconstructions 

(HI MART algorithm has been reconstructing constant fields with minimum error, but 
suffers smaller range of relaxation parameter (can give limited features on 
reconstruction). Anderson ART has spanned a greater range of relaxation parameter and 
given reconstructions with an error slightly more than GH MART. Error levels in Gilbert 
and Gordon ART are equal but more than Anderson and GH MART. Lent MART is 
again reconstructing constant fields with same range of relaxation parameter as GH 
MART but the error levels arc higher. Simple ART possibly can have slow as well as fast 
changes during reconstructions (Spanning whole possible range of relaxation parameter). 
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GBH MART and Maximum entropy method are found to have reconstructions with 

maximum error. 


Overall it can be said that GH MART is the best algorithm for constant fields followed by 
Lent MAR 1 . Minimum energy and Simple ART. Minimum energy method gives 
minimum error in step and impulse field reconstructions; whereas the error levels are 
high in GB1 1 MART and Maximum entropy method. 


Only one value of relaxation parameter is only used to reconstruct three fields using 
maximum entropy method. These algorithms are slow and reconstruction takes too long 
time. fable-(4. 22-4.24) shows the various reconstruction errors in maximum entropy 
method for three fields, reconstruction errors are more than minimum energy method. 


Constant field 


I 

/f. 

h H 

E c 

E d 

E, 

Iterations 

0 1 

i 2001 -01 

9.936E-02 

6.958E-0I 

1.G26E-0) 

4.099E-02 

11 


Titbie-422 Different error levels, number of iterations, and relaxation parameter are given for maximum 

entropy constant function reconstruction 


Step field 


T 

E a 

p 

■B 

E c 


Ee 

Iterations 



9.553E-02 

1.541E-01 

2.787E-02 

3.844E-02 

27 


T*bfe-4.Z£ Different error levels, number of iterations, and relaxation parameter are given for maximum 

entropy step function reconstruction 
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IabIe-4.24 Different error levels, number of iterations, and relaxation parameter are given for maximum 

entropy impulse function reconstruction 
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4.2 Simulated results 


0 


0.5 


i 


area average 0.5 



area average 0.5001 




Fig-4.9 Constant original field 


Fig-4.10 Reconstructed constant field 
using Simple ART 

(61 rays, 12 views, X = 1 .99) 


o 


05055 


I 


area average 05000 



05411 




area average 05002 




Fig-4. 1 1 Reconstructed constant field 
using Anderson ART 
<«1 rays, 12 views. k = 0.9 ) 


Fig-4.12 Reconstructed constant field 
using Gilbert ART 

(6 1 rays, 12 views, k = 0. 1 5 ) 
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area average 0.5002 



area average 0.4992 




Fig-4.13 Reconstructed constant field 
using Cordon ART 
(61 rays, 12 views, >1 = 0.15) 


Fig-4.14 Reconstructed constant field 
using GBH MART 
(61 rays, 12 views, k = 0.01 ) 



Flg-445 Reconstructed constant field 
using GH MART 
(61 rays* 12 views, k = 0.10 ) 


Fig-446 Reconstructed constant field 
using Lent MART 
(61 rays, 12 views, k = 0.22 ) 
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area average 03001 area average 03002 


Fig-4.|7 Reconstructed constant field 
sing Maximum entropy method 
<61 rays, 12 views, X = 0. 1 ) 


Fig-4. 18 Reconstructed constant field 
using Minimum eneigy method 
(61 rays, 12 views, X = 0.8) 



area average 0.725 area average 0.7251 


Fig-4, 19 Step original field Fig-4.20 Reconstructed step field 

using Simple ART 

(61 rays, 12 views, /l = 1.99) 
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Fig-4,21 Reconstructed step field 
using Anderson ART 
(61 ray*, 12 view*. X ~ 0.9) 


Fig- 4.23 Reconstructed step field 
using Gordon ART 
(61 rays, 12 views, X = 0. 1 5 ) 


0 



area average 0.7251 


Fig-4.22 Reconstructed step field 
using Gilbert ART 
(61 rays, 12 views, A = 0.15) 



0.8507 


M 



area average 0.7252 


Fig-4.24 Reconstructed step field 
using GBH MART 
(61 rays, 12 views, X = 0,0 1 ) 


0 



*iSXs: 


Fig-4.28 Reconstructed step field 
using Minimum energy 
(61 rays, 12 views, X = 0.8) 


Fig-4,27 Reconstructed step field 
using Maximum entropy 
(61 rays* 12 views, X = 0. 1 ) 


0.8587 


Hg-4.25 Reconstructed step field 
using GH MART 


(61 ravs, 12 views, A "0.1) 


area average 0.7248 


area average 0.7251 


Fig-4.26 Reconstructed step field 
using Lent MART 
(61 rays, 12 views, X = 0.22 ) 


area average 0.7251 


area average 0.7250 


0.9063 


0 




alea average 2.5044 


Mg-4.31 Rcconitrtictcd Impulse field 
using Anderson ART 
(61 rays, 12 views, X = 0.9 ) 


0 


m 


jg|g 

3.1086* 



area average 2.5040 


Fig-4.30 Reconstructed impulse field 
using Simple ART 
(61 rays, 12 views, X — l . 99 ) 



area average 25037 


Fig-4.32 Reconstructed impulse field 
using Gilbert ART 

(61 rays, 12 views, X = 0. 1 5 ) 



area average 2.5037 area average 2.5092 


Fig-4.33 Reconstructed impulse field Fig-4.34 Reconstructed impulse field 

using Gordon ART using GBH MART 

(61 rays, 12 views, A = 0 . 15 ) (61 rays, 12 views, A = 0 . 01 ) 



area average 2.5033 area average 23032 


Fig-4.35 Reconstructed Impulse field 
using GH MART 
<61 rays. 12 views, /l = 0.1) 


Fig-4.36 Reconstructed impulse field 
using Lent MART 
(61 rays, 12 views, A — 0.22 ) 
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0 


3.0391 




area average 2.5012 




3.1080 


area average 2.5044 


Fifc-4* J7 Reconstructed Impulse field 
using Maximum entropy 
61 rays* 12 views, A ~ 0. 1 ) 


Fig-4.38 Reconstructed impulse field 
using Minimum energy 
(61 rays, 12 views, A = 0.8) 
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4.3 Discussions over simulated results 

All the mentioned algorithms have been extensively studied for the range of relaxation 
parameter and five aforementioned error norms for constant, step, and impulse 
reconstructions. First two of the Above-mentioned reconstructed fields, can be used to 
construct any function, in addition the third function, which is studied in the present 
work, is impulse. The basic purpose of considering this function in the present study is to 
see the sensitivity of these algorithms towards a pulse change in the field. 

As earlier mentioned, relaxation parameter acts as a filter function , means if relaxation 
parameter is small it signifies slow changes during reconstruction. A study of relaxation 
puumctcr with RMS error with step function reconstructions is presented. It can be seen 
easily in the given I ; ig-(4. 1-4.8) that RMS error is minimum for fast changes, moreover it 
depends upon the type of field one want to reconstruct., except in GBH MART, GH 
MARI', and Minimum energy method, where medium changes are preferred. Anderson 
AR f and Simple ART have almost covered the whole possible range of relaxation 
parameter. Simple ART always shows a convergence whereas Anderson ART starts 
showing diverging characteristics as relaxation parameter goes beyond 1.1, and after this 
ii is at the verge of divergence. Gilbert ART and Gordon ART can be used up to 
relaxation parameter 0.15, beyond which they show divergence characteristics (i.e. not 
able to show the features corresponding to higher relaxation). GBH MART shows 
convergence up to 0.02 however this range goes a little above for both GH and LENT 
MARI' (up to 0.22). 

All the algorithms are highly sensitive to the initial guess, especially multiplicative ART 
methods. Initial guess can entirely change the error norms and their distributions for large 
relaxation p.uamefeis and may also lead to divergence. 

At last, it can be concluded that Anderson ART, GH MART, Lent MART and Minimum 
energy methods are the best algorithms for low error reconstructions. Simple ART and 
Anderson ART are good enough (can reconstruct maximum features of an image) for 
large range of relaxation parameter. GBH MART and Maximum Entropy method have 

shown larger error in reconstruction. 
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Chapter 5 


Two-fluid model for void fraction patterns in liquid metal two- 
phase flows 

I he uccui ate prediction oi two-phase flow dynamics with respect to the boiling and 
condensation conditions is essential for the safety of nuclear reactors. However for gas 
driven ADS sv stems two-phase flow studies, without interfacial mass transport are 

required 1 10}. 

1 tansient two-dimensional two-fluid model has been used to address the development of 
gas-driven ADS systems. I his model considers each phase separately expressed in terms 
of two sets of conservation equations governing the balance of mass, momentum and 
energy in each phase. I hey are complemented by the set of interfacial jump conditions 
(balance over discontinuities) with proper assumptions. Various suitable constitutive 
relations have been utilized (2. 8j. 

5.1 Model formulation and assumptions 

Model foinadniiou is now being discussed in detail [2, 8], This is a six-equation model, 
eon.'.idcrine each phase separately. Continuity equation has been written for both the 
phases separate!) with “smear densities” (the product of density and volume fraction). 
Interfacial mass transfer is absent in the present case (there is no phase change). 
Momentum equations are coupled with microscopic viscous momentum transfer, pressure 
equilibrium and interfacial momentum transfer. To simulate interfacial momentum 
transfer between two -phases in direct contact, a term “Momentum Transfer Coefficient 
( K )" has been considered and it is assumed that momentum transfer between the two- 
phases is proportional to the velocity difference between two-phases [2, 4]. Coupling 
interfaeial drag in this way has an advantage that it gives linearized implicit coupling 
between fields and K itself can vary with velocity [4]. 

Generalized interfaeial drag = K^Lh - L^| 

K ineorpoiatc.s steady and transient forces, in steady forces only skin drag has been 
Considered as turbulence has been neglected. On the other hand transient forces can be 
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fragmented in two parts, virtual mass force and Basset force. Virtual mass force is the 
force caused hv the virtual mass around the bubble. Basset force is caused by the 
acceleration of boundary' layer around the bubble. Faxen force has been neglected from 
our analysis because of its unimportance in bubbly, slug and chum turbulent flows. This 

needs special attention only in separated flows [12]. 

Certain elements of complex microphysics of exchange processes, such as the effect of 
momentum mixing (turbulence), have been neglected. Thus form drag has also been 
dropped out of steady drag formulations; moreover its associated energy dissipation and 
detailed distinction of interfacial properties from the bulk properties are also neglected. 
The effect of drag dissipation has been assigned completely to the heating of vapor [4]. 
interfacial forced convection heat transfer term, again for linear implicitness in numerical 
procedure, is written as [4], 

Interfacial forced convection term= R\(Ti - r x )| 

I he equations are also having microscopic viscous dissipation and microscopic 
conduction heat transfer within single phase [2, 4, and 8]. Turbulent heat flux and 
interfacial viscous heat dissipation (Faxen work) have been neglected [12], Due to the 
absence of interfacial mass transport only two jump conditions are required to handle 
interfacial discontinuities. Firstly the sum total of interfacial momentum is equal to zero 
and secondly, net interfacial heat transfer is equal to zero [8]. 

hirst jump condition gives rise to an assumption, which is the equality of mean pressures 
of the two phuws (i.e. /’ P„) [2, 4, and 8], However this is not a very good 

assumption, since with this assumption pressure waves travel with same speed in both the 
phases, i.e. gas and liquid. The isentropic acoustic velocities of liquid and gas differ a lot. 
This assumption has been frequently used earlier and hence has been utilized here too 
[17J. Second jump condition for interfacial heat transport is handled considering zero 
entropy production at interfaces with infinitely small interfacial thickness; 
implementation needs a higher value of R and leads to the equality of gas and liquid 
temperatures (2, 4]. 
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5.2 Conservation equations 

Continuity equations 


^-+v*te)=o 

Momentum equations 

'• + V' • (p s il s U J = = V/> + k(u , - U ,)+ v.(a £ j 

-</,)+ ?• (0 ■ -«k,) 


(5.1) 

(5-2) 


(5.3) 

(5.4) 


Note that the factors a and 1 -a (in equations 5.3 and 5.4) lie outside the pressure 
gradient. If' these factors were within the gradient, as some authors have proposed, then 
the equations could not represent static equilibrium for an inhomogeneous distribution in 

the absence of gravity [4], 

Hncrgy equations 




v4’/r.) 


• + y • te))+ 4i, -&} + 4; - 

- /| ‘2:?). + V . ((1 -a\l, )j + r{i\ ■ •/;)+ (1 - ak, • (W/ )+ v • A',(l - a^T, 


(5.5) 

(5.6) 


5.3 Constitutive equations 


Few more equations are needed to complete this model. These equations are as follows. 

(5J) 

=( 1 -#)# ( 5 . 3 ) 

Where, /> f ',/) k , are macroscopic liquid and gas densities (mass per unit of mixture 
volume) and are microscopic liquid and gas densities (mass of liquid per unit of 

liquid volume). Liquid is assumed incompressible and ideal equation of state tor gas has 

been utilized. 
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liquation of state for gas 


P = 


P*RT. 


m, 


(5.9) 


foi ideal gas under adiabatic reversible compression, i.e. isentropic process, acoustic 

velocity for gas is given as. 


a' 


tIP JP.tlV P „ V 2 y y 1 

. * ) = L PV = L nRT = r J~ RT 

dp c/i tip r m tn m mw 


a R = JyRT /m w = JyPlp 

Where, is the molecular weight of gas, and m is the mass of the gas. 


(5.10) 


The viscous stress tensors ay and a K are the usual Newtonian ones for liquid and gas 

respectively, in two-dimensional cylindrical coordinates, viscous stress tensors for liquid 
is being given and a similar set for gas holds [13J. 


{(T lir ) = Ip, 0*'- + X, 


i d(rui) dv, 
r dr g z 


W,r:)=Ml 


dv, du , 


dr dz) 

[a,..) -2(i, % + X, 


i d(ru,)d v L 

( r dr dz ) 


(<j, m ) = 2(i, u /,+X, 


i d(nt,),dvj 

l r dr dz j 


(5.11) 

(5.12) 

(5.13) 

(5.14) 


(t, is the liquid coefficient of viscosity (dynamic viscosity), and X, is the second 
coefficient of viscosity for liquids. The two coefficients of viscosity are related to the 

coefficient of hulk viscosity k by the expression [13], 


k = -(i,+X, 


(5.15) 


In general, it is believed that k is negligible except in the study of the structure of shock 
waves and in the absorption and attenuation of acoustic waves. For this reason, bulk 
viscosity k is ignored. With* = 0 , the second coefficient of viscosity becomes, 
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(5.16) 


4 


2 

3 


Mi 


Same expressions of second coefficient of viscosity have been utilized with gas also. 


5.4 Numerical procedure 

The basic numerical procedure is an extension of the Implicit, Continuous-fluid, Eulerian 
technique (ICE), which allows for fluid dynamic studies at all flow speeds [3]. Although 
the technique has been applied to generalized Eulerian-Lagrangian representations, the 
present multiphase flow version is restricted to a purely Eulerian mesh of computational 
cells. This technique is referred as an Implicit, Multi-Field (IMF) solution method [4], 
Among its properties are the following [4], 


1. Implicit treatment of mass convection and equation of state, so that flow speeds can 
range from far subsonic (incompressible) to supersonic. 


2. Implicit coupling among the fields, to allow forces that range from very weak to those 
strong enough to completely tie the fields together. 


The Eulerian mesh of finite-difference cells is used. In cylindrical coordinates these cells 
are toroids about the axis, rectangular in cross-section and with dimensions# and & 
whereas in Cartesian coordinates these cells are rectangular. 


Sr 



Flg-5.1 Layout of field variables and indices for a computational cell 

jwrr TTT'f-TrT 
HP StJffTXji 
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5.4.1 Convective fluxes 

I here are various ways in which mass, energy, and momentum convective fluxes can be 
written, the precise term being irrelevant to the conceptual formulation of IMF technique, 
but nevertheless crucial to the assurance of accuracy and numerical stability of the 
calculations. Partial donor cell ditlerencing has been utilized in this analysis, according to 

it (3|: 
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(a. + A /„) 
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<0 (5.16) 


(5.17) 


In the same way energy convection fluxes have been discretized. Partial donor cell flux is 
not space centered and has low order truncation errors. These contribute a positive 
diffusive effect and hence tend to automatically stabilize the numerical calculations. 

5.4,2 Finite difference formulations 

The finite difference formulations for all the equations have been presented below [4], 


(aJ'"'./ =(/’k)"u ~4(/A'v)' HI 'V -(p« u x r T''-l') lr > Sr 

-0>>,)"V ;)/& (5 - 18) 


*s,k ‘ (5 - I9) 
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(5.20) 


U'*v g y",.„ i = (G) -Sta n ",.,+L(p n +','^ - p n * ] , ,)/ Sz 

i «** utv , l -(v s r',, /t i] 


it' 1 ) '■< = (p')'.i -<5if((/7/«,r)' ,+l , + I / - (yo'/K,r) ,+l ,-i ; )//*,<5r 

( s - 21 ) 

(/’ «(£),♦;., -4- a " +l 4^" +, ' +1 '' ~p n+, ^) /Sr 

u'itK", t i } ,,{fa g '} 'it-',/ -(»/) ,,+1 '+^/) (5.22) 


rV,.*, «(/.) , 

» /,M n 2 


r\,)isz 


(5.23) 


Where, the quantities (O*),, 1 ,(<r) , , (L) t J , and (l) ' account for the effects of 

* M 4 , 2* ; ' 2 

momentum .. .-mcem-m gravity, and viscous stress. These have been used at time level n 
[4 j making this technique semi-implicit. 


(<;).. « ip * u K - ( p ' s " 2 /)' 1 - 1 ) ! r (+ \ 5r 

2 * 2 

-(p« u x v J‘*'r>-' 2 ) /Sz 

+ )'m., - r,or\, (cr^ )/ r , £r 

2 

-Sta n ,PJa m y^lJr , 


(5.24) 
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Where. 
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Whore. 


0’4, -(^rr)M n+I '^'-< + '^) /Sr 

+ {<*;,. ),.,(m/" +/ m (.,4/ + + u"*‘ ,-Lj-l )/ 45z 

+(*J, f (v; , ^w., t 5+v;^4-v^U,4 + vr^-/,-|)/^ 


( 5 . 33 ) 



The viscous stress components are discretized using standard central differencing. 
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1 he discrepancies in continuity equations can be written as [4], 


kk*k r )"-j J 

(k*k -(0>«kTj 

7 i 


1 


St / 8z 


St I r, Sr 


(5.43) 




(5.44) 


l here is no separate equation for pressure calculations; so pressure has been adjusted 
iteratively using Wwton-Rapshon method over continuity equations in order to keep the 
values of (/> I and (D i ) i , below a specified convergence criterion in each cell [3, 4 


and 7J, If void-fraction is more than 0.1 ( a ( ), pressure has been iterated over liquid 
pressure correction equation otherwise gas pressure correction equation is used in 
iterativ e explicit way, i.e. implicit. The value of a 1 does not affect the nature of solution 
although it may have some effects on convergence rate. Convergence will be faster when 
it will he iterated over liquid continuity equation, because liquid has higher acoustic 
velocity. The final purpose of pressure iteration is to obtain a pressure field for which the 
d: .v icn.tn.'k". in continuity equations are at minimum (compatible with convergence 
criterion). However instead of achieving full convergence in single iteration, (l - fi) 

times of previous iteration values of (D,\ ; and are achieved. In the present work 

over-relaxation is used {/I = 1.5). 

( \>n\ ei pence criterion is given below. 




w 'U'l, 

for al«x v ' 

< 

lO’ip'X 

for a"j > a ( 


(5.45) 
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Niwton kupshon Equation or liquid and gas pressure correction equations for pressure 

correction are given as. 


=(p,y -(Dr 


or 


(>'J" ={rj u 


a(A),„ 

l 

v ap ^ , 


(5.46) 


Now. it' accomplish the derivation of 


a( ^U and r a (4>, 


dP 

\ '•/ J 


V dP u J 


it is considerably more 


coincident to differentiate a slightly different expression for (D,)"*‘and(z) R )j’ +l , namely 

the foim obtained when the mass convection is derived from a strict centered 
differencing, with no propagation of donor-cell fluxing [4]. The converged solution is 


independent of 




and 


f'lA)„ 


■* 


dr. 


, it is entirely acceptable to use any desired set 


'■< J 


of saines for this coefficient, provided that the resulting convergence rate is not adversely 
affected, In many cases, it is sufficient to calculate the values only at the start of the 
iteration, and keep them fixed throughout the cycle at those initial values. However this 

arbitrariness does not apply to the specification of (D , )"*' and {D k )'*' , in this case, forms 
must be exactly as specified in Hq-(5.43) and (5.44). Now using Eqs-(5.19), (5.20), 


(5.22), (5.23), (5,4.3), and (5.44) 
below |4{, 


^U and f^l 


dP 


/./ j 


\ dP u j 


are calculated and are given 
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Where. 




J / sent topic 


(5.49) 


U,l, =( 1 /( 1 -<')>? 


(5.50) 


The two above mentioned equations are the general forms for calculating 

W !/>,1 ',“(/>.) > 

: „ ( and' 1 , ■ . However a simplified form of above equations has been used 

in the proem analysis as these forms affect the convergence and not the solution. 

When void traction is less than the critical value then iterations are carried out over liquid 
equation- t iquid metals possess very large values of (A , ) ; ; and (AT),,. The simplified 

' (/’ ) 

fe:m ot : can be written as [4L 

■ P/* : 

/ / . . . \ 

(5.51) 



f / 

' (pi)', Sr 2 Sz- " 

l j 

= // 

2$t[Sr 2 + &')^ 


1 ! ’ ■ ■ iik' change in velocity due to change in pressure in gas, from Eq-(5.47), the 

expression fur i is calculated as [4], 

j 


«'(/>.. ) ) , t 

" 1/U) 


r ,a , + r i« i 

V i 2 1 ' 2 2 y 


(<») 2 /r, (<&•)’ 


j (if , + « I 


w /(&)’ 




(5.52) 


Where, 


dp 

f 

V v K /f ' f J imitropic 

U,),. ,*( y <‘h 2 


(-'«), (sgj- 


(5.53) 

(5.54) 
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io calculate the void fraction in each cell, if a" ] <a c then 



(5.55) 

Andifer ’' ' a then 


<;'=i -fet'teX, 

(5.56) 

1° handle the interfacial energy transfer implicitly intermediate temperatures are defined 

like* 


f,=v ■■+(/, -(/,)" )/(c,y 

(5.57) 

r ‘> t ; +(/,-(/, f'Mcd 

(5.58) 


I he interfacial heat transfer is equally divided between the calculations of 7 , 1, and 

final calculation of (/ t . , (/ J; . 

1 and / are calculated during explicit iterative pressure calculations only considering 
enenp caused by interfacial drag and half of interfacial heat transfer. With known I g and 
/ , densities and temperature of both the phases can be calculated using their equation of 

states. 

flu* terms for /, and l, are given as. 
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( KDRAG) = A',”, 



(160) 


(/. ) . . - 2 (<■,)■ 


(T,)l -('Xj '{C,X, 

-b.XAhXAC'X 


M, W" - s <k, 

* « ' );, ),, /(c, ); /fac,^ (a E + «*j, ) 


/2 




(5.61) 


Final solution for internal energies is obtained, which accounts for the effect of 
convection. viscous work, single fluid conduction, and the remaining half of the 
interfaeia! heat transfer. Final expressions for internal energies are given below. 


(U:;' =(<;): 


[/, +i 2 + / 3 ] 


(5.62) 


Where. 


( (/’, t ' f f ((a- f - fa t ) 

(rtf/^Xt'TTOO-tCGO+tCGj) 
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(5.64) 
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(«;.)= -W,> & 

2 ' 2 

(co,)=fcX.“‘" Ml-fc U /& 
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(«r«,-2).(kr^U,-Ur^U 

v ~ 2 

+ c((w , “r v L;, -((«r»r , '-Li,) / '; 

- 4X-; 44' 4 ■ 44, j 


/4 


(5.65) 

(5.66) 


(5.67) 


(r<;,)=- r , (AT,)r , k ) itu -(r,)>.* (5S8) 


'V' <V, 


(<r, J -- r , (a:, )" , ; ( ((?„ ) fJ - (r, > 'A 


A/e 




(4,4;l+(UA4, 


(5.69) 


( 2 (C,)',(aI* , ) /2 < 5 - 70 > 


4 


fo-iCftL (w:/-«J 

.(A/&X(CT/ , 0'3)-(CG s ) + (CGj) fo*- ) 
(A/AX(CW4)-(CG 7 ) + (CG 8 )) 

+ (A ! r ,% ~ a ”j )i,j ( u i )i j ; 


(5.71) 
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( 5 . 79 ) 


-I Pd.,) 

- (St / & X(CI-7»G3) - (CG 5 ) + (CG b )) - (St / Sr\(C VPG 4) - (CG 7 ) + (CG 8 )) 

* X! 


V, t,’ - (< '/ - **;, |r, - (r, ), j + (r, ), , /(c, > - (J, ),, i(c, X, | 

(-(<':)" (P:T! +a K,) 

< ((< : )" *K, (i, t! iff, V(2(c, f" (a t' + #K, ) (S-W) 


The quantities (rr„ ) and (<t,) tire evaluated using time level (n+1) velocities and the time 

level (n) viscosities. 
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5.4.3 Solution procedure 


The solution procedure step-by-step is given as [4], 

1 . First set up the initial conditions for P , T g , T,, a , u g , u n v g , and v, in each cell 

of computational domain and some constant variables. 

2. Calculate initial values of p g , p,, p g , p t , I g , I,, C g , and C, for each cell in 

computational domain from gas and liquid equation of state (Eqs-(5.96)-(5.102), 

l''qs-(5.7)-(5.8)). 

3. Assign all above evaluated values to n th time level values, and calculate n n time 
level A\. , K l , p g , p, , and acoustic velocities for gas and liquid for each cell. 


4. Apply velocity and pressure boundary conditions, like no-slip, free-slip, constant 
pressure, specified-influx, continuous outflow. 

5. Calculate!/,), (/<),,, C (G),^,, and (G),,, + I using n th time level quantities, 

and calculate n‘ h time level K and R (Eqs-(5.24)-(5.27)). 

6. Solve the four momentum equations explicitly to calculate the eight velocities, 
and then using these velocities and n th time level densities calculate mass fluxes 

slqs-t5.39H5.42)). 

7. For each time step update, from /" to 

8. For each cell, till all cells have \(D,) or|(/)J<£ simultaneously. 


H, I* f 


P(D,)) (*<£«] 

'■v )' nr , 


for each cell using n' h time level quantities (Eqs-(5.51)- 


(5 52 )} 

HJ. t \ilculnh' (D, ) far a", < . and(D g )for <, > a2 (Eq-(5.44)). 

ft 3, If |(/>, I or |(/), ] < * go to Mfiral < OC C , and to 8.8 forCC^ > • 

$4 If j|/) ; | c)r [^ | > e calculate new pressure using Newton-Rapshon method with over- 

n'hvuttitm (relaxation parameter- 1.5) (Eq-(5.46)). 
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8.5. 


( 'alcidat t 


’ (h) and ( 7 /) 


omitting temporarily the effects of convection, work, and single-phase 


1 1 induct ton. Evaluate p, and 7) from liquid equation of state, and p g and T g from gas 

equation of state (Eqs-(5.59) and (5.61)). 

8.6. Again calculate the eight velocities namely, (u/) n+ j , (u/)” + j , (v, )" +l , , (v / )* +1 |> (u , 

l+ 2 J l ~2' J ’ S ' + v' 

y i ■ ( v , ; tlml (v s r i f rom momentum equations (Eqs-(5.39)-(5.42)). 

i ’ 1 * «. s 


& calculate mass fluxes for gas. 

8. ”, /. Si dve the gas continuity equation for [p g )"*' (Eqs-(5. 18)). 

8.-.:. Calculate al r, 

8 . "J < -adulate (/> )"'/ = (i -a"' \p, ), y 
iV. " /" !(/) ) go to the next cell 

& **„ 1 /f ( / ) , | > /; calculate new mass /luxes using the new [p t j and velocities and return to 8.1 

; vr r^vu'/vvvwc P f *' estimate 

S.S. b u " '> a calculate mass fluxes for liquid. 


H t H. l S> d \ « * //**» thpttil continuity equation far (/?, / (Eqs-(5> 2 1)), 

8.8.:, i die, date a", : : 1 ” (a )", /(/>/ • 

#•«.* tabulate (p'J/ =-<</(p?\, 

8.8.4 It :(/> ). r: go la the next cell. 

8. HJ, If j > /; calculate new mass fluxes using the new [p g )'/ and velocities and return to 8.1 

tar new pressure P ** 1 estimate* 

9. Calculate ( /J Hl and (l,)"'! implicitly considering the effects of convection, 

fa * j ^ ^ 1 / 

work, and single-phase conduction, and half of interfacial heat transfer. Evaluate 
Pl and from liquid equation of state, and p g and (t^ from gas 

equation of state (Eqs-(5.62) and (5.80)). 
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10. Calculations for one time step update are over return to 7 for next time step, and 

repeat same. 

5.4.4 Interfacial transfer formulations 

Interfacial drag coefficient and heat transfer coefficient are the two important parameters 
governing the actual interface transfer processes and need special attentions. The 
formulation for interfacial drag coefficient in undistorted particle regime (Viscous 

regime l has been given below [8, 11, and 12]. 

, v a K , 4. 2°. \MsL f d JL 

K V n J. ■' VCD (5.81) 

Where, /•'. , /> . /•• and fi m are the standard drag force, volume of a typical particle, 

virtual mass force, bubble radius and mixture viscosity, respectively [8, 11, and 12]. In 
the above formulation first term deals with steady force whereas last two terms deal with 
transient forces. Basically, the first term consists of interfacial skin and form drag. 

1 hmev cr the second term is virtual mass force, when a bubble accelerates then the virtual 
liquid mass around it also accelerates with the same velocity. Virtual mass force itself 
consists of two parts, one is transient and other is due to convection of liquid around the 
bubble Convection term for virtual mass force has been neglected because of high liquid 
viscosity. Third term is Basset force, which arises from acceleration ol viscous boundary 
layer around the bubble. Boundary layer formation around the bubble has also been 
neglected; therefore Basset force has been dropped out. An additional force, which is 
possible, ts Faxon force (a viscous contribution), which is the correction applied to virtual 
mas?, effect and Basset force to account for fluid velocity gradient. Faxen force has not 
been ■ in the present formulation. 

Hie bubbly flow regime is expected to exist for low void fractions, nearly less than 0.3. 
Assuming spherical bubble of uniform size the number of bubbles per unit volume is, 

n h =3a/(4m- h J ) ^ 5 ‘ 82 ^ 


r h , is the radius of a typical bubble. 
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(5.83) 


Intertacial area per unit volume in bubbly flow is defined as [8], 

A , = 3a / r h 

Drag resistance of a single vapor bubble can be written as 


=~ C oPA 


V\ 


(5.84) 


\\ here ( is the resistance coefficient, p is the liquid density, A is the typical projected 

area of a bubble and V is the relative velocity vector. 

t hi further simplification. 


A- A ) 


(5.85) 


So. the resistance or friction force per unit volume is obtained by multiplying Eq-(5.85) 

In ». |!2|, 

. . f’ 

{v,~y e ) (5.86) 
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• »■ -- ~-C n p, V, - V 

n, * 


Drat* resistance coefficient for the spheres is defined as [12J, 

0.42 


... 24 [l'f {).15Re a687 ) + 


* Re 
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1 4- 4.25 x 1 0 4 Re L16 


(5.87) 


he Rc\ n«»U!s number is given as. 


Re 


V 7 / 


K-K 


Pi 


(5.88) 


I::!.-:: h.i.e drag coefficient A', caused by steady forces, is defined as the interphase 
faction value divided by the difference of phase velocities. 



(5.90) 
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tS r. 


P; 


f)-V x 


24 f. 


\'imial mass force is given as [12J, 


aF ± 


^ I M a Pl 


d k (v k ) _ 

■ ! -v t( .vv, 


Dt 


(5.91) 


' • "• (* . ~ * $ ) ’ s relative velocity between phases, and C VM virtual mass coefficient. 

Where, the substantial derivative is given as, 

/iLLiJ+j? ,v( ) 

Dt dt K w 

tin further *t:t :j *Ii fie.iti. -n, by neglecting the convective part of virtual mass force and 

considering purely transient part, term becomes. 


ah\ 


C ni ap, 


dt 


(5.92) 


\’i:tual mass coefficient for bubbly How is given as [8], 


C 


l\l 



l + 2a 

1 -a 


1 hen h.tc: pit t .e virtual mass coefficient A', is defined as. 
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■ a 

2 

**** 

1 ■ f 2a 

1-a . 

ap, 

dt 


K- 




K = K,+K 2 


(5.93) 


(5.94) 

(5.95) 


20 

Presently very large value of interlacial heat transfer coefficient (A =1.0x10 ) has 

been taken. It has been assumed that interface does not control the interfacial heat 
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been found quite efficient for bubbly flows to handle temperature discontinuities. 

5.4.5 Numerical stability 

Stewart (W>) considered the numerical stabilities of two-fluid model in a quite 
dhusified manner [5J. It was said that two-fluid modeling of two-phase flow may yield a 
svskm ot partial differential equations having complex characteristics; this results in a 
mathematically ill-posed initial value problem. The behavior of finite difference 
equations with and without exchange terms has also been explored by Stewart [5]. Two- 
lluid finite difference formulations without interfacial exchange term have been 
considered first and it was concluded that if the convective limit (courant number [6], 
with connection propagation, not sonic propagation) on the time step is obeyed, the 
highest frequency modes {wavelength small multiple of cell dimension) should not grow 
;v**mo!:iv.i!l\. even for the ill-posed two-fluid problem. Now turning to low frequency 
mode*, time and space step should not be too small. 

Vv itti momentum exchange considerations, it was shown that the product of momentum 
exchange coefficient and space step should not be too small. Moreover, partial donor cell 
differencing of convective fluxes contributes mitigating positive diffusive term and 

makes the numerical technique stable enough. 

At present in bubbly flow, having large momentum exchange, suitable value of space 
step has been taken to assure stability. Appropriate time step has been taken (neither too 
large nor too small 1 based upon convective courant number because of explicit handling 
of convective terms in momentum equations [6, 5). 
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5.5 Modeling of LBE-Argon system 

liquid metal two-phase flow system considered is a lead-bismuth and argon gas two- 
phase flow. Calculations have been carried out in a two-dimensional Cartesian coordinate 
system, dividing the ~QxJ) square cm test field into 23x23 grids, including fictitious 
cells, Ml the diseiiti/ed equations governing the complex microphysics of two-phase 
flow hav e been solved using point-by-point gauss siedel method with over relaxation in 
piessuie eot tee lion using New ton-Rapshon method. LBE (Lead-Bismuth Eutectic) and 
aigon gas are entering at the bottom and are leaving out from the top of the physical 
domain, a square obstacle has been considered at the centre of the physical domain. Its 
Schematic has been given below (Fig-5.2). 


Schematic of the physical domain 

Physical and computational Domain 

luleriun cells have been used with cell sizes AX = 0.95 cm, AY = 0.95 cm and 
\i 1.5 * 10 ’sec. Time step has been taken based upon iterative pressure corrections. 
In ensuring that with a single iteration during pressure correction, fluid should not travel 
more than one cell (courant stability criterion). Iterative explicit technique during 
pressure collection facilitates nearly implicit treatment of momentum equation for sonic 
pi. >p, motion However, convective terms are handled explicitly and therefore limit the 
stability on convective propagation. Thus ICE allows higher time step than pure explicit 
technique (3 1. Calculations have been done without gravity. 
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5.5.1. Boundary conditions 

I he boundary conditions have been enforced by considering fictitious cells surrounding 

the computing mesh. 

5.5.2. No-slip boundary condition on a wall 

For rigid walls, the tangential velocities in fictitious cells are set by the reflection with a 
change in sign. However, normal velocity has been taken zero. 

5.5.3. Continuous outflow boundary 

Continuous out- flow boundaries are gradient free in their normal direction. 

5.5.4. .Specified inflow boundary 

X and Y components of velocities are defined. 


For initial conditions, the whole physical domain has been considered mostly filled with 
I, BF (umi fraction 0.05) at l.Ox l0 6 dyne$/cm 2 pressure and 403 K temperature. At t=0 
sec liquid and gas both are introduced at the bottom with v=10 cm/sec, u=0.0 cm/sec at 
1 ,{) Itf d\nes enf and 403 K temperature with void fractions=0.3. 

Due to undistorted particle or viscous regime consideration, bubble radius has been 
assumed constant throughout the flow and value used in all calculations is 0.20 cm. For 
simplicity, the interfacial heat transfer coefficient R is set to large constant values i.e. 
ver> tight coupling is assumed in terms of heat transfer, (i.e. temperatures of both liquid 
and gas tire equal ) or in other words entropy production at the interface is assumed zero 
p, 4). High values of R and K work ideally with highly dispersed flows, but if the flow 
regime changes from dispersed to separate or some other regime then this is not a good 
assumption. In that case, both R and K need special attentions. 


5.6 Thermo-physical properties of the LBE-Argon system 

ReUimi..hips for the Lead-Bismuth physical properties are given below [47], 

Density (gtriU-m*)- 10.737 - 1 .375 x 10 3 T JM 


Itermai 


? i * 

4 D* el * i < 


{ergs / sec cm 0 C ) — 7.26 x 1 0 +1. 23x10 T IBl , : 


(5.96) 

(5.97) 
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Specific heat (ergs 1 gm ” c) = 1 .465 x 1 0 6 


(5-98) 


Dynamic viscosity (gw/cmsec) = 3.26xl0" 2 -6.26x\Q~ 5 T [BE +4.63x10 ~ S T^ BE (5.99) 
Relationships lor the Argon gas physical properties are given below [48, 49], 

Density [gm an ')~ (P Rl)m w (5.100) 

Phermal conductivity [ergs / see an" c) 

(1032.74 t 4.4747‘„ - 1.725 x 10“ , 7’; /< + 5.528 x lO" 7 ^ ) (5.101) 

Specific heat [ergs t gm" (’)=■■ 5.207x1 0' 1 (5.102) 


Ihe relationship of dynamic viscosity of Argon gas with temperature has been utilized 

from PERRY'S CHEMICAL ENGINEER’S HANDBOOK [49]. 

Reference Tempcmture (Kehins) 403 
Reference pressure oh nes cm' 1 / x / if 

Where, temperatures T, H or T w are in "C , and pressure (P) in dynes/cm 2 . 
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Because of utilizing staggered grids in the formulations, both the velocity components are 
not defined at the same location. So to plot velocity vectors linear interpolations have 
been done to define both the velocity components at centre of the cells. Fig-6.5-6.8 
represents the velocity vectors of both the fields at three time levels. 
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Fig-6.6 Gas velocity vectors at t=1.5 sec 
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Fig-6.8 Gas velocity vectors at t=4.5 sec 


|. ig.,, am,. I ; represents the liquid velocity vectors in the physical domain at three time 
levels. I torn all above presented velocity vectors it car. be seen that enough equality exist 
between the gas ami liquid velocities, thus verifies the actual dynamics and enough 

realistic formulations ol imcrl'acial drag of bubbly (lows. 
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Both the fluids are entering at the bottom and show smooth velocity profiles, as time 
i no t eases the iomiation of eddies in the corner locations takes place and velocity vectors 
are compatible enough with void fraction patterns for realistic bubbly flows. 


20 


20 


E 

u 


•? 

$ 

>- 


0 


0 




X-axis (cm) 20 



>- 



^ 1 ' * , , r f 



X-axis (cm) 20 


Fig»63> Liquid velocity vectors at t=0.75 sec 


Fig-6.10 Liquid velocity vectors at t=1.5 sec 
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Fig-6. 1 1 Liquid velocity vectors at t -3.0 sec 


Fig-6. 12 Liquid velocity vectors at t=4.5 sec 


Solution of energy equations, with assumption of no entropy production at interfaces also 
has agreement with actual bubbly flows. Small temperature variations down the flow 
have been reported. Temperature values of both the fields are nearly same and are in 
agreement with bubbly Hows with higher R. 

Initial conditions play a crucial role in actual two-phase flow stabilities; these may be 
uniform or non-uniform but should not have discontinuities. LBE is a high-density fluid 
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with strong hydrostatic pressure variations, so for tall structural simulations in liquid 
metal specifically, non-uniform initial conditions are recommended. Presently uniform 
initial conditions have been employed in absence of gravity. 
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Chapter 7 


Conclusions and future scope 

Gamma ray attenuation technique with suitable tomographic inversion algorithms and 
two-fluid simulations of liquid metal two-phase flows are having different result sets. 
Two-fluid modeling has shown superiority in terms of result sets; along with void 
distributions it gives pressure and velocity profiles. Conclusions on both the techniques 
based upon the present work have been drawn. 

7.1 Tomographic algorithms 

First conclusions are drawn over tomographic algorithms suitable for corrupted void 
fraction projection data, based upon the studies carried out it can be said that GH MART, 
Lent MART, minimum energy and Simple ART algorithms have given encouraging 
results for constant field reconstructions. 

Minimum energy method gives minimum error in step and impulse field reconstructions, 
and spans enough range of relaxation parameter (almost whole possible range). The error 
levels in GBH MART and maximum entropy method have been reported at maximum. 
Initial guess of the field values in these tomographic algorithms have their own 
importance; different initial guesses may give different error levels in reconstructions. 
Studies related to initial guess are due and can be extended in future work. 

Gamma-ray attenuation technique can be used to collect projection data on ADS leg 
followed by their inversion in void fraction distributions. 

7.2 Two-fluid simulations 

In addition to void fraction distributions, two-phase flow modeling gives velocities, j 
pressure drops, etc. Two-fluid modeling of LBE-Argon in a simple geometry has enough 
agreement with realistic pure bubbly flows. Still the experimental verification is required 
to verify the results of two-fluid simulations for its further application over gas-driven ; 

ADS target systems. ; 
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In the present literature, formulations of interfacial transport terms are just limited to pure 
bubbly flow; this limitation can be extended to apply the same model with same solution 
procedure over other flow regimes (like slug, churn turbulent, annular flows), with 
change in interfacial coupling terms. This is a two-fluid single pressure model; and can 
be updated to a two-pressure model to deal gas and liquid separately in terms of pressure. 
Initial conditions also play a crucial role in real two-phase stabilities; non-uniform initial 
conditions should be employed without any discontinuity. For tall structural simulations, 
non-uniform initial conditions are desirable because of strong hydrodynamic pressure in 
liquid metals. No coalescence or disintegration of the bubbles has been accounted. The 
bubble tracking equation can also be coupled with present six field equations to account 
for coalescence etc. In interfacial heat transfer calculations, tight coupling between two 
fields, in agreement of pure bubbly flows, is assumed with high interfacial heat transfer 
coefficient (very small interfacial thickness, no entropy production at interfaces). But if 
flow regime changes from bubbly to annular etc., then proper analytical model for 
interfacial heat transfer should be considered. Turbulence modeling can also be added in 
the present formulation to extend the model for turbulent two-phase flows. 

Here implementations have been done on a simple geometry for demonstration purpose; 
however studies are due on actual ADS target module geometry considering gravity and 
transient inflow boundary conditions to accommodate the actual gravity dominating ADS 
loop flow. 

At last, it can be said that gamma ray attenuation technique with suitable tomographic 
inversion algorithms can also be used to verify the results of two-fluid simulations. 
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APPENDIX A 


Implementation of optimization reconstruction algorithms 

A.l Minimization of energy through Lagrangian multiplier optimization technique 
1 he / projection, which is a path integral, is approximated as 

N 

i = l, ,M (A.1.1) 

/ = ! 

Minimum energy functional is defined as 




(A. 1.2) 


Now based upon Lagrangian multiplier optimization technique, 


M 

F =-L/, ! -La 

/=! 


i=i 


I lore, X, are hi unknown Lagrangian multipliers. 
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Substituting Hq-(A. 1 .4) in Eq-(A. 1.1) and simplifying 


(A. 1.3) 
(A. 1.4) 


/sAf 

zL/ ^ kL ~ ~^ P k (A. 1 .5) 

i=! ;=l 

Eq-(A.1.5) is a system of linear equations with the multipliers X, as unknown. This set of 
equations can be solver through Gauss-Seidel method for X, , and thereafter using A,. , field 
can be reconstructed. 
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A.2 Maximization of entropy through Lagrangian multiplier optimization technique 
Maximum entropy functional is defined as 


J 


(A.2.1) 


Now based upon Lagrangian multiplier optimization technique for optimization of a 
function under some constraint, 




;=1 


lL M ’ufr p , 


7=1 


(A.2.2) 


1 lore. A t are M unknown Lagrangian multipliers. 


On simplification. 



r / jsM 

f, = exp 1 - 


W, 


0/ 


/ = 1 


Now putting /, from Hq-(A.2.4), in Eq-(A.l.l), 


£ FI ( ex p(- X ' w 'j )) =p > ex p(0 

/ I /al 


Replacing exp(~/l, ) = x, in Eq-(A.2.5), 


Assuming, 






^„m^= p t ^) 
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(A.2. 3) 


(A.2.4) 


(A.2. 5) 


(A.2. 6) 


(A.2. 7) 


&,(*;)= ^ ex P0) 


(A.2. 8) 


90 



Now Eq-(A.2.6), can be solved by Gauss-Seidel method after linearization using 
truncated Taylor series expansion, 

Let x'" be the first guess for the solution set x ( (n) . 


/’ exp( 1) = £, (x, )= g, (x', <n )- 


K dx i A(» 


(xP ) -xf ) ) + 


\ dx lJ 
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( 4 2) -^) 


+. 




fit'' 




(A.2.9) 


Now using Bq-(A.2.9) x, 1 " 1 can be calculated, and then Lagrangian multipliers to 

calculate field values. 
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